/*
 * Example_CalciteSolubility.h
 *
 *  Created on: 26 Oct 2011
 *      Author: allan
 */

#ifndef EXAMPLE_CALCITESOLUBILITY_H_
#define EXAMPLE_CALCITESOLUBILITY_H_

// C++ includes
#include <iostream>
#include <iomanip>
using namespace std;

// GeoReact includes
#include <GeoReact.h>

void Example_CalciteSolubility()
{
	GeoManager manager;
	
	manager.SpeciateAqueousSpecies({"H2O", "Na[+]", "Cl[-]", "Calcite"});
	manager.SetMineralSpecies({"Calcite"});
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	cout << system << endl;
	
	GeoEquilibrator equilibrator(system);
	
	equilibrator.ImposeChargeBalance();
	equilibrator.ImposeSpeciesMass("H2O", 55.508);
	equilibrator.ImposeElementMass("Na[+]", 2.0);
	equilibrator.ImposeElementMass("Cl[-]", 2.0);
	equilibrator.ImposeSpeciesMass("CO2(g)", 1.13610543);
	equilibrator.ImposeSpeciesMass("Ca[2+]", 1.0E-8);
	equilibrator.ImposeSpeciesMass("Calcite", 1.99820162);
	
	GeoState state(system);
	
	state.SetTemperature(60.0);
	state.SetPressure(80.0);
	
	equilibrator.Equilibrate(state);
	
	cout << state << endl;
	
	GeoReactor reactor(system);
	
	const double hour = 3600;
	const double day  = 24 * hour;
	
	double t = 0.0, tEnd = 30*day;
	
	cout << setprecision(10) << fixed;

	while(true)
	{
		cout << setw(20) << left << t/hour;
		cout << setw(20) << left << state.Get("Moles", "Calcite")    * 100.09;
		cout << setw(20) << left << state.Get("Moles", "CO2(g)")     * 44.01;
		cout << setw(20) << left << state.Get("Moles", "CO2(aq)")    * 44.01;
		cout << setw(20) << left << state.Get("Moles", "CaCO3(aq)")  * 100.09;
		cout << setw(20) << left << state.Get("Moles", "NaHCO3(aq)") * 84.01;
		cout << setw(20) << left << state.Get("Moles", "NaCO3[-]")   * 83.0;
		cout << setw(20) << left << state.Get("Moles", "CaHCO3[+]")  * 101.09;
		cout << setw(20) << left << state.Get("Moles", "HCO3[-]")    * 61.02;
		cout << setw(20) << left << state.Get("pH") << endl;
		
		if(t == tEnd) break;
		
		reactor.React(t, tEnd, state);
	}
}

#endif /* EXAMPLE_CALCITESOLUBILITY_H_ */
